#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""
@FileName: brent.py
@Version: 1.0
@Date: 2023/6/21 17:28
@Description:
"""

from typing import Tuple
import math
from root_paras import Paras
from find_section import find_section

Tri = Tuple[float, float, float]
Quad = Tuple[float, float, float, float]


# tri: xi, a_, b 三元组
def sort_points(fun, tri: Tri):
    xi, a, b = tri
    ya, yb, yi = math.fabs(fun(a)), math.fabs(fun(b)), math.fabs(fun(xi))
    if not (yi < ya and yi < yb):
        xi: float = (a + b) / 2.0
        if fun(a) * fun(xi) < 0:
            b, xi = xi, (a + xi) / 2.0
        else:
            a, xi = xi, (xi + b) / 2.0
    return xi, a, b


def after(fun, quad: Quad) -> Tuple[bool, Tri]:
    xi, a, b, c = quad
    a_new, b_new = a, b
    if fun(a) * fun(c) < 0:
        b_new = c
    else:
        a_new = c

    # 后向误差 f(xi) 得到改进 and 包含根的区间减少一半
    succ = False
    if math.fabs(fun(c)) <= math.fabs(fun(xi)):
        xi = c
        sect = b - a
        if b_new - a <= 0.5 * sect or b - a_new <= 0.5 * sect:
            succ = True
        if b_new - a < sect:
            b = b_new
        elif b - a_new < sect:
            a = a_new
    xi, a, b = sort_points(fun, (xi, a, b))
    return succ, (xi, a, b)


def ioi(fun, tri: Tri, sol_tol):
    xi, a, b = tri
    q = fun(a) / fun(b)
    r = fun(xi) / fun(b)
    s = fun(xi) / fun(a)
    numer = r * (r - q) * (xi - b) + (1.0 - r) * s * (xi - a)
    denom = (q - 1) * (r - 1) * (s - 1)
    term2 = -numer / denom

    c = xi + term2
    return after(fun, (xi, a, b, c))


def regul_falsi(fun, tri: Tri, sol_tol):
    xi, a, b = tri
    dy: float = fun(b) - fun(a)
    if math.fabs(dy) < sol_tol:
        print(f'Secant method, dy = {dy} ~> 0')

    c = (a * fun(b) - b * fun(a)) / dy
    return after(fun, (xi, a, b, c))


def bisect(fun, tri: Tri, sol_tol):
    xi, a, b = tri
    c = float(a + b) / 2.0
    return after(fun, (xi, a, b, c))


def brent(fun, paras: Paras = Paras()):
    a, b = find_section(fun, paras)
    if ret := fun(a) * fun(b) > 0:
        raise RuntimeError(f'invalid section [{a}, {b}], fun(a)*fun(b)={ret} > 0!')
    if b < a:
        b, a = a, b

    sol_tol = paras.sol_tol
    sec_tol = paras.sec_tol
    xi: float = (a + b) / 2.0

    # lucky!
    if math.fabs(fun(xi) - 0) <= sol_tol:
        return xi, (a, b)
    if math.fabs(fun(a) - 0) <= sol_tol:
        return a, (a, b)
    if math.fabs(fun(b) - 0) <= sol_tol:
        return b, (a, b)

    if b - a < 2.0 * sec_tol:
        b = a + paras.dft_sect
    xi, a, b = sort_points(fun, (xi, a, b))

    n: int = 0
    loops = paras.loops
    while b - a > sec_tol:
        n += 1
        if n > loops:
            raise RuntimeError(f'Secant iteration: loop over 1000 times, with interval {a - b}')

        # lucky!
        if math.fabs(fun(xi) - 0) <= sol_tol:
            break

        succ, (xi, a, b) = ioi(fun, (xi, a, b), sol_tol)
        if not succ:
            succ, (xi, a, b) = regul_falsi(fun, (xi, a, b), sol_tol)
            if not succ:
                succ, (xi, a, b) = bisect(fun, (xi, a, b), sol_tol)

    if math.fabs((val := fun(xi)) - 0) <= sol_tol:
        print(n)
        return xi, a, b
    else:
        raise RuntimeError(f'Secant iteration: did not converge to zero, instead -> {val}')


if __name__ == '__main__':
    fun = lambda x: -x ** 3 + 80 * x - 1
    res = brent(fun)
    print(res)
